/**
 * 
 */
package exactDiagonalisation;

import org.apache.commons.math.stat.StatUtils;

import Jama.EigenvalueDecomposition;
import Jama.Matrix;

import common.Stat;
import common.Tuple;

/**
 * @author vladimir
 */
public class Hamiltonian {
	int N;
	double minE = 0;
	int nStates = -1;
	int nSites = -1;
	int dimension = -1;
	double J = -10000000;
	int[][] neighbors = null;
	double[] Eig_norm = null;

	/* function takes number of state and returns coordinates on the lattice */
	/**
	 * private int[] i2xy(int i) { int[] result = new int[2]; result[0] = i %
	 * this.N; result[1] = (int) (i / this.N); return result; }
	 **/
	/* function takes number of state and returns state as array of +1 and -1 */
	public int[] n2state(int i, int tN) {
		char[] binstr = Integer.toBinaryString(i).toCharArray();
		int[] tresult = new int[tN];
		for (int j = 0; j < binstr.length; j++) {
			tresult[j + tN - binstr.length] = Integer.parseInt(String
					.valueOf(binstr[j]));
		}

		int[] result = new int[tN];

		for (int j = 0; j < tN; j++) {
			result[j] = 2 * tresult[j] - 1;
		}

		return result;
	}

	/*
	 * function takes state and returns state number *
	 */
	public int state2n(int[] state) {
		int tN = state.length;
		int[] tresult = new int[tN];

		for (int j = 0; j < tN; j++) {
			tresult[j] = (int) ((state[j] + 1) / 2);
		}
		StringBuilder stresult = new StringBuilder();
		for (int j = 0; j < tN; j++) {
			stresult.append(Integer.toString(tresult[j]));
		}

		return Integer.parseInt(stresult.toString(), 2);

	}

	public int xy2i(int x, int y, int N) {
		return y * N + x;
	}

	private int sum_neighbors(int site, int[] state) {

		int result = 0;
		for (int i = 0; i < this.i_neighbors[site].length; i++) {
			result += state[this.i_neighbors[site][i]];
		}
		return result;
	}

	private Tuple[][][] xy_neighbors = null;
	private Tuple[][] x_neighbors = null;
	private int[][] i_neighbors = null;
	Matrix H;
	private EigenvalueDecomposition Eig;

	// private DenseMatrix H;

	/** Constructor for test purposes **/
	public Hamiltonian() {

	}

	/** Another Constructor for test purposes **/
	public Hamiltonian(int N) {
		this.N = N;
	}

	public Hamiltonian(int N, double J, double B, int[][] neighbors) {
		this.N = N;
		this.dimension = neighbors[0].length;

		this.nSites = (int) Math.pow(this.N, this.dimension);

		this.xy_neighbors = new Tuple[N][N][neighbors.length];
		this.x_neighbors = new Tuple[N][neighbors.length];
		this.J = J;
		this.neighbors = neighbors;

		// now lets build 2d array representing matrix
		this.nStates = (int) Math.pow(2, this.nSites);// number of
														// different
														// states
		double tLattice[][] = new double[nStates][nStates];

		if (dimension == 2) {
			for (int x = 0; x < this.N; x++) {
				for (int y = 0; y < this.N; y++) {
					for (int tn = 0; tn < neighbors.length; tn++) {
						int tx = x + neighbors[tn][0];
						if (tx == -1) {
							tx = this.N - 1;
						} else if (tx == this.N) {
							tx = 0;
						}
						int ty = y + neighbors[tn][1];
						if (ty == -1) {
							ty = this.N - 1;
						} else if (ty == this.N) {
							ty = 0;
						}
						this.xy_neighbors[x][y][tn] = new Tuple(0, tx, ty);
					}
				}
			}
			this.i_neighbors = new int[this.nSites][neighbors.length];
			for (int x = 0; x < this.N; x++) {
				for (int y = 0; y < this.N; y++) {
					for (int tn = 0; tn < neighbors.length; tn++) {
						this.i_neighbors[xy2i(x, y, this.N)][tn] = xy2i(
								xy_neighbors[x][y][tn].x,
								xy_neighbors[x][y][tn].y, this.N);

					}
				}
			}

			for (int i = 0; i < nStates; i++) {
				int[] state = n2state(i, this.nSites);
				double H_J = 0;
				for (int l = 0; l < state.length; l++) {
					int site = state[l];
					H_J -= site * sum_neighbors(l, state);
				}
				tLattice[i][i] = J * H_J / 2;
				for (int k = 0; k < this.nSites; k++) {
					state[k] *= -1;
					int j = state2n(state);
					state[k] *= -1;
					tLattice[i][j] = B;
				}
			}
		} else if (dimension == 1) {
			for (int x = 0; x < this.N; x++) {
				for (int tn = 0; tn < neighbors.length; tn++) {
					int tx = x + neighbors[tn][0];
					if (tx == -1) {
						tx = this.N - 1;
					} else if (tx == this.N) {
						tx = 0;
					}
					this.x_neighbors[x][tn] = new Tuple(0, tx, 0);
				}
			}
			this.i_neighbors = new int[this.nSites][neighbors.length];

			for (int x = 0; x < this.N; x++) {
				for (int tn = 0; tn < neighbors.length; tn++) {
					this.i_neighbors[x][tn] = this.x_neighbors[x][tn].x;
				}
			}

			for (int i = 0; i < this.nStates; i++) {
				int[] state = n2state(i, this.nSites);
				double H_J = 0;
				for (int l = 0; l < state.length; l++) {
					int site = state[l];
					H_J -= site * sum_neighbors(l, state);
				}
				tLattice[i][i] = J * H_J / 2;
				for (int k = 0; k < this.nSites; k++) {
					state[k] *= -1;
					int j = state2n(state);
					state[k] *= -1;
					tLattice[i][j] = B;
				}
			}
		}

		this.H = new Matrix(tLattice);
		this.Eig = new EigenvalueDecomposition(this.H);
		// double[] eigenvalues = Eig.getRealEigenvalues();

		this.minE = StatUtils.min(Eigenvalues());// minimum Energy
		this.Eig_norm = new double[this.nStates];// Energy specturum with
													// substracted minimum
													// energy

		for (int i = 0; i < Eigenvalues().length; i++) {
			this.Eig_norm[i] = Eigenvalues()[i] - this.minE;
		}
	}

	// returns double array with eigenvalues
	public double[] Eigenvalues() {
		return this.Eig.getRealEigenvalues();
	}

	/*
	 * as parameters eats, T - temperature and returns partition function
	 * multiplied exp(minE/T) by
	 */
	public double Z(double T) {
		double result = 0;
		// double minE = Stat.min(Eigenvalues());// minimum Energy
		for (double Energy : this.Eig_norm) {
			double power = -(Energy) / T;
			result += Math.exp(power);
		}
		return result;
	}

	/*
	 * eats temperature as parameter and returns Energy at this temperature
	 */

	public double Energy(double T) {
		double result = 0;
		// double minE = Stat.min(Eigenvalues());// minimum Energy
		for (double E : Eigenvalues()) {
			double power = -(E - this.minE) / T;
			result += E * Math.exp(power);
		}
		return result / Z(T) / this.nSites;
	}

	/* eats temperature and returns entropy per site at this temperature */
	public double Entropy(double T) {
		return Math.log(Z(T)) / this.nSites - this.minE / T / this.nSites
				+ Energy(T) / T;
	}

	/* eats temperature and returns specific heat per temperature */
	public double C(double T) {
		double E2 = 0;
		for (double E : this.Eig_norm) {
			E2 += E * E * Math.exp(-E / T);
		}
		E2 /= (Z(T) * this.nSites);

		return Math.pow(T, -2) * (E2 - Math.pow(Energy(T) * N, 2));
	}

	/*
	 * takes temperature as parameter and returns Magnetisation
	 */
	public double M(double T) {
		double result = 0;
		double[] M = new double[this.nStates];
		Matrix V = this.Eig.getV();
		for (int n = 0; n < this.nStates; n++) {
			M[n] = 0;
			for (int j = 0; j < V.getRowDimension(); j++) {
				double Sz = Stat.sum(this.n2state(j, this.nSites))
						/ this.nSites;
				M[n] += Math.pow(V.get(j, n), 2) * Sz;

			}
			result += M[n] * Math.exp(-this.Eig_norm[n] / T);
		}
		return result / Z(T);
	}

	/*
	 * takes temperature as parameter and returns Magnetisation squared
	 */
	public double M2(double T) {
		double result = 0;
		double[] M2 = new double[this.nStates];
		Matrix V = this.Eig.getV();
		for (int n = 0; n < this.nStates; n++) {
			M2[n] = 0;
			for (int j = 0; j < V.getRowDimension(); j++) {
				double Sz = Stat.sum(this.n2state(j, this.nSites))
						/ this.nSites;
				M2[n] += Math.pow(V.get(j, n), 2) * Math.pow(Sz, 2);
			}
			result += M2[n] * Math.exp(-this.Eig_norm[n] / T);
		}
		return result / Z(T);
	}

}
